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ABSTRACT 

Computations of cosmic-ray transport 
based upon finite-difference methods 
are afflicted by instabilities, inac- 
curacies, and artifacts. To avoid 
these problems, we have developed a 
Monte Carlo formulation which is 
closely related not only to the finite- 
difference formulation, but also to the 
underlying physics of transport phenom- 
ena. Implementations of this approach 
are currently running on the Massively 
Parallel Processor at Goddard Space 
Flight Center, whose enormous computing 
power overwhelms the poor statistical 
accuracy that usually limits the use of 
stochastic methods. In a Monte Carlo 
simulation of rectilinear transport, 
the coherent and diffusive effects that 
appeared are in good quantitative 
agreement with both finite-difference 
and analytic calculations. 
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INTRODUCTION 

The diffusion idealization, which has 
been almost universally invoked in 
discussions of cosmic-ray transport, is 
easy to treat analytically. However, 
many observed phenomena give clear evi- 
dence for non-diffusive effects. One 
example is the so-called "scatter free" 
propagation of kilovolt solar electrons 
(Ref. 1), which is inconsistent with 
diffusion, but which can readily be 
interpreted in terms of a coherent mode 
of propagation. This mode is novel, 
but it is just a manifestation in a 
dynamic situation of non-diffusive 
effects similar to those considered in 


the steady-state by classical transport 
theory (Ref. 2). Although these 
effects have been described analyti- 
cally in References 3 and 4, the theory 
is very complicated. Consequently, 
there is a need for reliable numerical 
computations which bypass these com- 
plexities and yield concrete results 
suitable for comparison with observa- 
tions. This paper explores such 
computations based on the well-known 
Monte Carlo method and compares them to 
computations based upon finite- 
difference methods. To limit the 
number of parameters, and to focus the 
discussion on computational methods, 
the present discussion is limited to 
the case of rectilinear one-dimensional 
propagation of cosmic-rays along a 
uniform magnetic guiding field on which 
are superimposed random fields. This 
leaves out the effects of adiabatic 
focusing by non uniform guiding fields 
and of convective motion of the back- 
ground medium, which are very important 
in the interplanetary context, but it 
includes two other essential aspects of 
charged particle transport. These are 
a strong inhibition of transport per- 
pendicular to the guiding field and a 
pronounced anisotropy of the pitch- 
angle scattering by random fields. 
Finally, note that the magnetic fields 
are visualized as static, which means 
that the velocities of individual 
particles are constant, and that there 
is no interaction among particles in 
the extremely tenuous distribution of 
cosmic-rays. This situation differs 
from those considered by plasma 
physics, but it is closely analogous to 
those treated by classical transport 
theory. 
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TRANSPORT EQUATIONS 

Under the circumstances outlined above, 
particle transport is described by 


from that described by the finite- 
difference formulation, has a signifi- 
cant effect that will be discussed 
below. 
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in which f is particle density in phase 
space, y is the cosine of the pitch- 
angle measured with respect to the 
guiding field, and z is distance paral- 
lel to the guiding field. The variable 
s = Vt, where V is particle velocity, 
plays the role of a temporal parameter. 
The Fokker-Planck coefficient of pitch- 
angle scattering is given by 
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where 1 is the mean free path, and q is 
an index that measures the anisotropy 
of scattering (Ref. 5). 


In the discrete formulation, the con- 
tinuous variables are replaced by a 
three-dimensional grid whose spacings 
are Az, Ay, and As = VAt, and the 
derivatives appearing in equation (1) 
are replaced by their finite— difference 
analogs. These replacements lead to 
the difference equation 


M = P M+^ f MH,Z " f M,Z^ 

+ p m--j (f M,Z “ f M-l,Z* 
+ ^M.Z-l “ f M,z)‘ 


which gives the change in the distribu- 
tion function during a temporal incre- 
ment As, and which can readily be 
solved by standard numerical methods. 
In the Monte Carlo formulation, the 
random history of a large number of 
particles is followed under the assump- 
tion that the coefficients Pm^/ 
appearing in equation (3) can be inter-2 
preted as probabilities that y will 
change by ±Ay in each time step. In 
this formulation, the particles move a 
distance y AS in each step. This 
slight difference in the evolution of z 


In Figure 1, the transition probabili- 
ties are plotted as a function of y for 
two different assumptions about the 
anisotropy of scattering. Both curves 
refer to the same mean free path, A - 
480. The square symbols describe 
strongly anisotropic scattering (q = 
1.8) similar to that occurring in 
interplanetary space. Evidently, the 
scattering near y = 0, which is the 

boundary between forward ( y > 0) and 

backward (y < 0) hemispheres, is much 
weaker than that occurring within each 
hemisphere. Consequently, this 

boundary acts as a physically signifi- 
cant barrier which particles find 
difficult to penetrate. In contrast, 
the circular symbols, which refer to 
isotropic scattering (q = 1.0) similar 
to that occurring in molecular colli- 
sions or neutron diffusion, describe 
relatively weak and nearly uniform 
scattering with no feature at y = 0. 



Figure 1. Scattering probabilities for 
X = 480 plotted as a function of pitch- 
angle cosine. The square points for 
anisotropic scattering define a curve 
that is qualitatively different from 
the one defined for classical isotropic 
scattering by circular points. 
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THE ALGORITHM 

To implement the Monte Carlo scheme 
outlined above, each particle was 
assigned an integer distance and an 
index corresponding to p that lies 
between 0 and 7 . Because these 
parameters occupy only three bytes, 
there was plenty of storage for several 
parallel arrays of particles. Conse- 
quently, the results given below are 
based on 31 arrays which contained 31 * 
16384 = 507904 particles. Their fate 
was determined by a single parallel 
array of random integers (ranging from 
-32767 to +32767) that was rotated at 
each temporal step relative to the 
fixed arrays of particle data, and 
updated every 31 steps. To implement 
changes in the pitch-angle cosine, 7 
positive integers were chosen in such a 
way that the probabilities of the 
current random numbers being larger 
than this integer are those given by 
Figure 1. Then the angular index was 
incremented for those particles whose 
current random integer was positive and 
greater than the appropriate proba- 
bility integer, and decremented for 
those whose random integer was negative 
and less than the probability integer 
with its sign reversed. This approach 
satisfies the basic requirement that 
the probabilities of incrementing, 
decrementing and leaving unchanged the 
pitch-angle must add to unity. After 
the pitch— angles had been updated, each 
particle's distance was changed 
accordingly. When the desired number 
of temporal steps had been carried out, 
particles were binned into an array 
according to distance and pitch-angle. 

ISOTROPIC INJECTION 

Results obtained from the MPP for 
isotropic injection of particles 
uniformly distributed over pitch-angles 
are presented in Figure 2 as plots 
versus distance of the total number of 
particles in each of 50 bins. This sum 
over pitch-angles is a measure of the 
isotropic particle density. Because 
the total number of particles was 



Figure 2. Density profiles shortly 
after an isotropic injection. The 
profile for isotropic scattering is 
featureless, while that for anisotropic 

scattering exhibits two prominent 

coherent pulses. 

large, statistical errors are small 
and, consequently, are not shown 
explicitly. However, slight irregu- 
larities in some parts of the curves 

given an indication of their magnitude. 

In Figure 2, the darker curve refers to 
anisotropic scattering (square symbols 
in Fig. 1), while the lighter one 
refers to classical isotropic scat- 
tering (circular symbols). These 

density profiles describe a situation 
very shortly after injection when the 
particles have had time to move a 
maximum distance of only one mean free 
path. The former curve exhibits two 
peaks, which are completey absent from 
the latter one. Qualitatively, this 

manifestation of the coherent mode 

appears because equal numbers of 

particles become nearly uniformly 

distributed in each hemisphere, while 
very few particles penetrate from one 
hemisphere to the other through the 
region of weak scattering near p = 0. 
Particles that stay together in each 
hemisphere, move with nearly the same 
velocity parallel to the field, but 

statistical fluctuations of individual 
velocities give rise to a peak centered 
around the average displacement. Such 
features are designated as coherent 
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pulses. There are two of them in 
Figure 2 because the injection was 
symmetrical. 

THE COHERENT MODE 

A detailed analysis (Ref. 3) predicts 
that the coherent pulses discussed 
above form a moving Gaussian density 
profile given by 


(z-V*t) 


F = 


(4irD*t) 


1/2 


(4) 


Consequently, to bring out these 
aspects more clearly, this section will 
focus on highly collimated injection 
into a region of anisotropic scatter- 
ing. More specifically, it will deal 
with injection of particles in the 
forward hemisphere at a single value 
y = 0.6. Figure 3 shows a density 

profile obtained under these circum- 
stances, but with all other conditions 
the same as in Figure 2. Evidently, 
the profile is very sensitive to 
initial conditions, for only one 
coherent pulse appears, and, in place 
of the second pulse, it is accompanied 
by a broad feature that is designated 
as the diffusive wake. 


where V* is a characteristic velocity, 
which is close to half the particle 
velocity, and D* is the coefficient of 
dispersion, which describes the broad- 
ening that arises from statistical 
fluctuations of individual velocities. 

Although the isotropic injection of 
Figure 2 is the most natural choice of 
initial conditions, it obscures certain 
important aspects of the coherent mode. 



Figure 3. Density profile after a 
collimated injection, with all other 
conditions the same as in Figure 2. 
Evidently, profiles depend sensitively 
upon conditions at injection. The 
dotted curve, derived from a finite- 
difference calculation, is in good, but 
not perfect agreement with the Monte 
Carlo results. 


The dotted curve gives the result of a 
finite-difference computation based on 
Equation 3. The overall agreement is 
excellent, but there is a slight 
displacement of the coherent peaks, and 
the dotted peak is slightly broader 
than the solid one. The first of these 
discrepancies is a trivial artifact 
arising from imperfections in the 
binning, but the second arises from an 
important difference between the two 
computations. More specifically, the 
finite-difference implementation leads 
to a dispersive effect which can be 
described by 

= -|v( Az - -|as) , (5) 

and the total dispersion is a super- 
position of the physical effect 
described by D* and the numerical 
artifact described by By. In the 
present example, the condition that 
ensures that physics dominates , 

D* > D //> i s met, for D* = 4D^, but the 
numerical dispersion is large enough to 
account for the slight discrepancy 
between the peaks in Figure 3. To 
confirm this interpretation, a finite- 
difference calculation was performed 
with D* = , and the expected further 

broadening of the coherent pulse was 
seen. These details are significant, 
for they illustrate that the Monte 
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Carlo method, by virtue of the differ- 
ence mentioned above in the evolution 
of z, is relatively free of the arti- 
facts which plague finite-difference 
calculations. (See Ref. 6 for a 
discussion of published work that has 
been affected by such problems.) 
Numerical dispersion is an especially 
insidious artifact, for it leads to 
results that seem plausible, but are 
quantitatively in error. 

ANGULAR DISTRIBUTIONS 

To bring out additional features of the 
coherent mode. Figure 4 shows further 
results from the anisotropic injection 
that was discussed above. Here, parti- 
cle density has been averaged sepa- 
rately over the forward and backward 
hemispheres. A comparison of the two 
curves reveals that particles in the 
pulse are overwhelmingly collimated in 
the forward direction while those in 
the wake are predominantly collimated 
backward. This behavior suggests that 
the wake can be interpreted as parti- 
cles scattered out of the pulse which 
subsequently move coherently away from 
it in the backward direction. 



Figure 4. Densities in the forward and 
backward hemispheres under the same 
conditions as in Figure 3. The 
coherent pulse contains particles 
moving forward. It is accompanied by a 
broad wake of particles moving 
backward. 


Detailed angular distributions tabula- 
ted at four points indicated on Figure 
4, and normalized to the same total 
number of particles, are presented in 
Figure 5. The distribution at the peak 
of the pulse (Curve B, circular points) 
is essentially a mirror image of the 
one at the center of the wake (Curve C, 
square points). Both describe near 
isotropy in one hemisphere associated 
with nearly zero intensity in the 
other. This distribution, which is 
characteristic of the coherent mode, 
adds weight to the interpretation 
outlined above. 

On the fringes of the density profile, 
at points A and D, angular distribu- 
tions (triangular points) are highly 
collimated in the forward and backward 
directions, respectively. In the 
interplanetary context, the fringes 
would correspond to the the very first 
particles to arrive in a solar event, 
which are particularly difficult to 
describe theoretically, but which are 
very important in connection with 
accurate timing of events on the sun. 
On the MPP , Monte Carlo methods, can 
give a useful description of this phase 
of solar events, but statistical fluc- 
tuations limit the applicability of the 
method on less powerful computers. 



Figure 5. Normalized angular distribu- 
tions at locations indicated in Figure 
4. These distributions exhibit a 
remarkable mirror symmetry. 
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THE DIFFUSION LIMIT 

Eventually, scattering destroys the 
strong anisotropies of the coherent 
mode, and the featureless wake becomes 
dominant. In this regime, where the 
familiar theory of diffusion is 
applicable, the density profile is 
described by a Gaussian whose width is 
controlled by the coefficient of 
diffusion D = X V/3 (Ref. 7). Figure 6 
shows a profile computed by the MPP in 
this regime for Vt = 64 * X (solid 
curve). Evidently, this profile is in 
very good agreement with the dotted 
curve which is based upon diffusion 
theory. 



Figure 6. Density profiles in the 
diffusion regime. Monte Carlo results 
are in very good agreement with the 
dotted Gaussian derived from diffusion 
theory. 

CONCLUSIONS 

Results obtained from the MPP with the 
aid of Monte Carlo methods are equiva- 
lent in every detail to those based 
upon careful use of more traditional 
methods, but they are less subject to 
error and are closer to the physics. 
These characteristics offer tremendous 
advantages in the investigation of 
exotic transport regimes for which no 
theoretical description is available. 
In particular, the formulation of 
problems in which particles gain or 
lose energy leads to prohibitively 
large conventional computations, but 


their Monte Carlo versions are not sig- 
nificantly more complicated than the 
one described here. We intend to ex- 
ploit these advantages in the investi- 
gation of two such problems: Adiabatic 

deceleration of cosmic-rays due to 
expansion of the solar-wind, and the 
loss of energy by electrons in radio 
sources due to synchrotron radiation. 
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